#include "HY1C_out.h"
/* =========================================================== */
/* Module l2_hdf_generic.c                                     */
/*                                                             */
/* Functions to open and write a multi-sensor (generic) l2     */
/* file in HDF format.                                         */
/*                                                             */ 
/* Written By:                                                 */
/*     Bryan A. Franz, SAIC GSC, March 1998.                   */
/*     Gary Fu,        SAIC GSC, March 1999.                   */
/*     Joel M. Gales,  Futuretech, Sept. 1999.                 */
/*     Gene Eplee, SAIC GSC, SeaWiFS Project, December 2000.   */
/*           Update time correction and mirror side            */
/*                       factors.                              */
/*     Gene Eplee, SAIC, SeaWiFS Project, March 2004.          */
/*           Convert time correction and mirror side           */
/*           factors to simultaneous exponentials.             */
/* =========================================================== */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h> 
#include <time.h>
#include <math.h>
#include "l12_proto.h"
#include "l2_hdf_generic.h"
#include "hdf_utils.h"
#include "time_utils.h"
#include "mfhdf.h"
#include "l2prod.h"

#define ROUND(x) ((((x) >= 0)?0.5:-0.5) + (x))

/* Global variables to facilitate communication */
static float bad_float  = BAD_FLT;
static float fill_float = FIL_FLT;
static float fill_int   = FIL_INT;
static float fill_byte  = FIL_BYT;
static int32_t  Lambda    [NBANDS+NBANDSIR], *Lambda_p;
static int32 numScans;
static int32 numPixels;
static int32 numBands;
static int32 numBandsIR;
static int32 spix; 
static int32 cpix; 
static int32 epix;
static int32 cscan; 
static int32 nctl; 
static int32 *ictl;
static int32 *jctl;

static FILE *fp_meta = NULL;

/* -------------------------------------------------------- */
/* -------------------------------------------------------- */
#define RFACTOR 100
int32 get_ctl(int32_t ctl_pt_fact, int32 ictl[], int32 jctl[])
{
    int32 nctl;
    int32 dctl;
    int32 i;

    if (ctl_pt_fact <= 0) {

        if (numPixels < RFACTOR) {
            nctl = numPixels;
            for (i = 0; i < numPixels; i++) {
               ictl[i] = i+1;
            }

        } else {
            dctl = (int32_t) numPixels/RFACTOR;
            nctl = MIN((int32_t) numPixels/dctl + 1,numPixels);
            for (i = 0; i < (nctl-1); i++) {
               ictl[i] = i * dctl + 1; 
            }
            ictl[nctl-1] = numPixels;
         }

    } else {
        dctl = ctl_pt_fact;
        nctl = MIN((int32_t) numPixels/dctl + 1,numPixels);
        for (i = 0; i < (nctl-1); i++) {
           ictl[i] = i * dctl + 1; 
        }
        ictl[nctl-1] = numPixels;
    }

    for (i = 0; i < numScans; i++)
        jctl[i] = i+1;

    return(nctl);
}

/* -------------------------------------------------------- */
/* -------------------------------------------------------- */
int16 *float2int16(float32 fbuf[], int32_t spix, int32_t npix, int incr, 
                 float slope, float offset)
{
    static  int16 *ibuf = NULL;
    float32 fval;
    int32_t    i;

    double maxval = slope*  32767  + offset;
    double minval = slope*(-32767) + offset;

    if (ibuf == NULL) {
        if ( (ibuf = calloc(numPixels,sizeof(int16))) == NULL) {
            fHY1C_out(stderr,
            "-E- %s line %d: Unable to allocate buffer space.\n",
            __FILE__,__LINE__);
            exit(1);
        }
    }

    for (i=0; i<npix && i<numPixels; i++) {
        fval = fbuf[spix + i*incr];
        if (fval >= maxval)
            ibuf[i] =  32767;
        else if (fval <= minval)
            ibuf[i] = BAD_INT;
        else
            ibuf[i] = ROUND((fval - offset)/slope);
    }

    return(ibuf);
}
        

/* -------------------------------------------------------- */
/* -------------------------------------------------------- */
uint8 *float2uint8(float32 fbuf[], int32_t spix, int32_t npix, int incr, 
                 float slope, float offset)
{
    static  uint8 *ibuf = NULL;
    float32 fval;
    int32_t    i;

    double maxval = slope*255 + offset;
    double minval = slope*  0 + offset;

    if (ibuf == NULL) {
        if ( (ibuf = calloc(numPixels,sizeof(int8))) == NULL) {
            fHY1C_out(stderr,
            "-E- %s line %d: Unable to allocate buffer space.\n",
            __FILE__,__LINE__);
            exit(1);
        }
    }

    for (i=0; i<npix && i<numPixels; i++) {
        fval = fbuf[spix + i*incr];
        if (fval >= maxval)
            ibuf[i] = 255;
        else if (fval <= minval)
            ibuf[i] = BAD_BYT;
        else
            ibuf[i] = ROUND((fval - offset)/slope);
    }

    return(ibuf);
}
        

/* -------------------------------------------------------- */
/* -------------------------------------------------------- */
float *float2float(float32 fbuf[], int32_t spix, int32_t npix, int incr, 
                 float slope, float offset)
{
    static  float32 *fbuf2 = NULL;
    float32 fval;
    int32_t    i;

    if (fbuf2 == NULL) {
        if ( (fbuf2 = calloc(numPixels,sizeof(float32))) == NULL) {
            fHY1C_out(stderr,
            "-E- %s line %d: Unable to allocate buffer space.\n",
            __FILE__,__LINE__);
            exit(1);
        }
    }

    for (i=0; i<npix && i<numPixels; i++) {
        fval     = fbuf[spix + i*incr];
        fbuf2[i] = MAX((fval - offset)/slope,BAD_FLT);
    }

    return(fbuf2);
}
        

/* -------------------------------------------------------- */
/* -------------------------------------------------------- */
VOIDP scale_sds(float *data, l2prodstr *p, int nbands)
{
    VOIDP pbuf;
    int   band = MIN(MAX(p->prod_ix,0),nbands-1);
    int   npix = p->dim[1];

    switch (p->datatype) {
      case DFNT_INT16:
        pbuf = (VOIDP) float2int16(data,band,npix,nbands,p->slope,p->offset);
        break;
      case DFNT_UINT8:
        pbuf = (VOIDP) float2uint8(data,band,npix,nbands,p->slope,p->offset);
        break;
      case DFNT_FLOAT32:
        pbuf = (VOIDP) float2float(data,band,npix,nbands,p->slope,p->offset);
        break;
      default:
        fHY1C_out(stderr,
        "-W- %s Line %d: Unknown data type %d.\n",__FILE__,__LINE__,p->datatype);
        pbuf = data;
        break;
    }

    return (pbuf);
}


/* -------------------------------------------------------- */
/* -------------------------------------------------------- */
float *unscale_sds(VOIDP *data, l2prodstr *p, int32_t spix, int32_t npix, int incr)
{
    static float32 *fbuf = NULL;

    float32 fval, *fptr;
    int16   ival, *iptr;
    uint8   bval, *bptr;
    int32_t    i;

    if (fbuf == NULL) {
        if ( (fbuf = calloc(numPixels,sizeof(float32))) == NULL) {
            fHY1C_out(stderr,
            "-E- %s line %d: Unable to allocate buffer space.\n",
            __FILE__,__LINE__);
            exit(1);
        }
    }

    switch (p->datatype) {
      case DFNT_UINT8:
        bptr = (uint8 *) data;
        for (i=0; i<npix && i<numPixels; i++) {
	  bval = bptr[spix + i*incr];
          fbuf[i] = bval*p->slope+p->offset;
        }
        break;
      case DFNT_INT16:
        iptr = (int16 *) data;
        for (i=0; i<npix && i<numPixels; i++) {
	  ival = iptr[spix + i*incr];
          fbuf[i] = ival*p->slope+p->offset;
        }
        break;
      case DFNT_FLOAT32:
        fptr = (float32 *) data;
        for (i=0; i<npix && i<numPixels; i++) {
	  fval = fptr[spix + i*incr];
          fbuf[i] = fval*p->slope+p->offset;
        }
        break;
      default:
        fHY1C_out(stderr,"-W- %s Line %d: Unknown data type %d product %d.\n",
		__FILE__,__LINE__,p->datatype,p->cat_ix);
        break;
    }

    return (fbuf);
}


/* -------------------------------------------------------- */
/* Assign SDSes to Vgroups                                  */
/* -------------------------------------------------------- */
int MakeVgroups(filehandle *l2file){

  int32	h_id;
  int32 v_id;
  int32 sd_id = l2file->sd_id;
  int   i;

  /* Do we have the extra meta-data for SeaWiFS */
  int seawifs_meta = 0;
  if (l2file->sensorID == SEAWIFS) {
      int32 sds_id;
      if ( sd_select(sd_id,"scan_ell",&sds_id) == 0 )
          seawifs_meta = 1;
  }

  h_id = Hopen(l2file->name, DFACC_RDWR, 0);
  if(h_id == FAIL){
    fHY1C_out(stderr,"-E- %s line %d: Hopen() failed for file, %s.\n",
    __FILE__,__LINE__,l2file->name);
    return(HDF_FUNCTION_ERROR);
  }
  Vstart(h_id);

  /* Scan-Line Attributes */
  DPTB( v_attach(h_id, &v_id)			);
  Vsetclass(v_id, "Per File Data");
  Vsetname(v_id, "Sensor Band Parameters");
  DPTB( AddSdsToVgroup(sd_id, v_id, "wavelength")	);
  DPTB( AddSdsToVgroup(sd_id, v_id, "vcal_gain")	);
  DPTB( AddSdsToVgroup(sd_id, v_id, "vcal_offset")	);
  DPTB( AddSdsToVgroup(sd_id, v_id, "F0")	  	);
  DPTB( AddSdsToVgroup(sd_id, v_id, "k_oz")		);
  DPTB( AddSdsToVgroup(sd_id, v_id, "Tau_r")		);
  Vdetach(v_id);

  if (seawifs_meta) {

    /* Sensor Tilt */
    DPTB( v_attach(h_id, &v_id)			);
    Vsetclass(v_id, "Per File Data");
    Vsetname(v_id, "Sensor Tilt");
    DPTB( AddSdsToVgroup(sd_id, v_id, "ntilts")		);
    DPTB( AddSdsToVgroup(sd_id, v_id, "tilt_flags")	);
    DPTB( AddSdsToVgroup(sd_id, v_id, "tilt_ranges")	);
    Vdetach(v_id);

  }

  /* Scan-Line Attributes */
  DPTB( v_attach(h_id, &v_id)			);
  Vsetclass(v_id, "Per Scan Data");
  Vsetname(v_id, "Scan-Line Attributes");
  DPTB( AddSdsToVgroup(sd_id, v_id, "year")		);
  DPTB( AddSdsToVgroup(sd_id, v_id, "day")		);
  DPTB( AddSdsToVgroup(sd_id, v_id, "msec")		);
  DPTB( AddSdsToVgroup(sd_id, v_id, "slon")		);
  DPTB( AddSdsToVgroup(sd_id, v_id, "clon")		);
  DPTB( AddSdsToVgroup(sd_id, v_id, "elon")		);
  DPTB( AddSdsToVgroup(sd_id, v_id, "slat")		);
  DPTB( AddSdsToVgroup(sd_id, v_id, "clat")		);
  DPTB( AddSdsToVgroup(sd_id, v_id, "elat")		);
  DPTB( AddSdsToVgroup(sd_id, v_id, "csol_z")		);
  Vdetach(v_id);

  /* Geophysical Data */
  DPTB( v_attach(h_id, &v_id));
  Vsetclass(v_id, "Per Scan Data");
  Vsetname(v_id, "Geophysical Data");
  for (i=0; i<l2file->tot_prod; i++) {
    DPTB( AddSdsToVgroup(sd_id, v_id, l2file->l2_prod_names[i])	);
  }
  Vdetach(v_id);

  /* Navigation */
  DPTB( v_attach(h_id, &v_id)			);
  Vsetclass(v_id, "Per Scan Data");
  Vsetname(v_id, "Navigation Data");
  DPTB( AddSdsToVgroup(sd_id, v_id, "longitude")	);
  DPTB( AddSdsToVgroup(sd_id, v_id, "latitude")		);
  DPTB( AddSdsToVgroup(sd_id, v_id, "cntl_pt_cols")	);
  DPTB( AddSdsToVgroup(sd_id, v_id, "cntl_pt_rows")	);
  DPTB( AddSdsToVgroup(sd_id, v_id, "tilt")		);
  if (seawifs_meta) {
      DPTB( AddSdsToVgroup(sd_id, v_id, "orb_vec")	);
      DPTB( AddSdsToVgroup(sd_id, v_id, "sun_ref")	);
      DPTB( AddSdsToVgroup(sd_id, v_id, "att_ang")	);
      DPTB( AddSdsToVgroup(sd_id, v_id, "sen_mat")	);
      DPTB( AddSdsToVgroup(sd_id, v_id, "scan_ell")	);
      DPTB( AddSdsToVgroup(sd_id, v_id, "nflag")	);
  }
  Vdetach(v_id);


  Vend(h_id);

  if(Hclose(h_id) != SUCCEED){
    fHY1C_out(stderr,"-E- %s line %d: Hclose(%d) failed for file, %s .\n",
    __FILE__,__LINE__,h_id,l2file->name);
    return(HDF_FUNCTION_ERROR);
  }
  return(LIFE_IS_GOOD);
}


/*-------- -------------------------------------------------------- */
/* Open an HDF file and store global attributes in it.              */
/* ---------------------------------------------------------------- */
int openl2_hdf( filehandle *l2file)   
{
    int32_t   status;
    float  Gain  [NBANDS];
    float  Offset[NBANDS];
    float  Fonom [NBANDS];
    float  Fobar [NBANDS], *Fobar_p;
    float  Tau_r [NBANDS], *Tau_r_p;
    float  k_oz  [NBANDS], *k_oz_p;
    char   title[256];
    char   soft_id[200];     /* software version info */
    int32_t   sd_id;
    int32_t   sds_id;
    int32_t   tot_prod;
    int    i, k, ii;
    l2prodstr *p;
    char   *tmp_ptr, tmp_str[2048];
    int32_t    n;
    VOIDP   pbuf;
    float  *fbuf;
    char    wavestr[5];
    char    prefix[256];
    char    suffix[256];
    char   *p1, *p2;
    int32_t    nbvis;

#ifdef SEADAS
    char seadas_vs[64], data_center[128], prog_version[64];
    char station_name[128], station_code[5];
    float station_lat, station_lon;
#endif

    if (strcmp(l2file->input->metafile,"") != 0) {
        fp_meta = fopen(l2file->input->metafile,"w");
        if (fp_meta == NULL) {
            fHY1C_out(stderr,
            "-E- %s line %d: Unable to open specified meta-data file, %s .\n",
            __FILE__,__LINE__,l2file->input->metafile);
            return(HDF_FUNCTION_ERROR);
	}
    }   
	    

    /* Create the HDF file */
    sd_id = SDstart(l2file->name, DFACC_CREATE);
    if(sd_id == FAIL) {
       fHY1C_out(stderr,
       "-E- %s line %d: Could not create HDF file, %s .\n",
       __FILE__,__LINE__,l2file->name);
       return(HDF_FUNCTION_ERROR);
    }   
    l2file->sd_id = sd_id;
   
    /* Get number of bands, pixels, scans */
    numPixels = l2file->npix;
    numScans = l2file->nscan;
    numBands = l2file->nbands;
    numBandsIR = l2file->nbandsir;
    spix     = 0;
    cpix     = numPixels/2;
    epix     = numPixels-1;
    cscan    = numScans/2;

    /* Get control-point array */
    if ( (ictl = calloc(numPixels,sizeof(int32))) == NULL) {
        fHY1C_out(stderr,
        "-E- %s line %d: Unable to allocate control-point array.\n",
        __FILE__,__LINE__);
        return(MEMORY_ALLOCATION_ERROR);
    }
    if ( (jctl = calloc(numScans,sizeof(int32))) == NULL) {
        fHY1C_out(stderr,
        "-E- %s line %d: Unable to allocate %d scan control-point array.\n",
        __FILE__,__LINE__,numScans);
        return(MEMORY_ALLOCATION_ERROR);
    }
    nctl = get_ctl(l2file->ctl_pt_incr,ictl,jctl);
    
    /* Get sensor-specific attributes */
    if ((n = rdsensorinfo(l2file->sensorID,l2file->input->evalmask,NULL,NULL)) != numBands) {
        fHY1C_out(stderr,"-E- %s Line %d:  Error reading sensor table. %d %d\n",
        __FILE__,__LINE__,numBands,n);
        return(-1);
    }
    rdsensorinfo(l2file->sensorID,l2file->input->evalmask,"Lambda",(void **) &Lambda_p);
    rdsensorinfo(l2file->sensorID,l2file->input->evalmask,"Fobar", (void **) &Fobar_p    );
    rdsensorinfo(l2file->sensorID,l2file->input->evalmask,"Tau_r", (void **) &Tau_r_p );
    rdsensorinfo(l2file->sensorID,l2file->input->evalmask,"k_oz",  (void **) &k_oz_p  );

    for (i=0; i<numBands; i++) {
        Lambda [i] = Lambda_p[i];    
        Fobar  [i] = Fobar_p [i];   
        Tau_r  [i] = Tau_r_p [i];   
        k_oz   [i] = k_oz_p  [i];
        Gain   [i] = l2file->input->gain[i];
        Offset [i] = l2file->input->offset[i];
        if (l2file->input->outband_opt >= 2) 
            get_f0_thuillier_ext(Lambda[i],BANDW,&Fonom[i]);
        else
	  Fonom[i] = Fobar[i];
    }
    for (i=numBands; i<numBands+numBandsIR; i++)
        Lambda [i] = Lambda_p[i];    


    /*                                                                  */
    /* Build list of SDS names from input parameters                    */
    /* ---------------------------------------------------------------- */
    /*                                                                  */
    tot_prod = prodlist(l2file->sensorID,l2file->input->evalmask,l2file->l2prod,
                        l2file->def_l2prod,l2file->l2_prod_names);

    strcpy(l2file->l2_prod_names[tot_prod++],"l2_flags"); 

    HY1C_out("\n\nThe following products will be included in %s.\n",l2file->name);
    for (i=0; i<tot_prod; i++)
      HY1C_out("%d %s\n",i,l2file->l2_prod_names[i]);

    l2file->tot_prod = tot_prod;

    /*                                                                  */
    /* Create the scan-line SDSes                                       */
    /* ---------------------------------------------------------------- */
    /*                                                                  */
    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "year",                                          /* short name      */
    "Scan year",                                     /* long name       */
    NULL,                                            /* standard name   */
    "years",                                         /* units           */
    1996, 2038,                                      /* valid range     */
    0,0,                                             /* slope           */
    DFNT_INT32,                                      /* HDF number type */
    1,                                               /* rank            */
    numScans, 1, 1,                                  /* dimension sizes */
    "Number of Scan Lines", NULL, NULL               /* dimension names */
    ) );

    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "day",                                           /* short name      */
    "Scan day of year",                              /* long name       */
    NULL,                                            /* standard name   */
     "days",                                         /* units           */
    1,366,                                           /* valid range     */
    0,0,                                             /* slope           */
    DFNT_INT32,                                      /* HDF number type */
    1,                                               /* rank            */
    numScans, 1, 1,                                  /* dimension sizes */
    "Number of Scan Lines", NULL, NULL               /* dimension names */
    ) );

    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "msec",                                          /* short name      */
    "Scan-line time, milliseconds of day",           /* long name       */
    NULL,                                            /* standard name   */
    "milliseconds",                                  /* units           */
    0,86399999,                                      /* valid range     */
    0,0,                                             /* slope           */
    DFNT_INT32,                                      /* HDF number type */
    1,                                               /* rank            */
    numScans, 1, 1,                                  /* dimension sizes */
    "Number of Scan Lines", NULL, NULL               /* dimension names */
    ) );

    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "tilt",                                          /* short name      */
    "Tilt angle for scan line",                      /* long name       */
    NULL,                                            /* standard name   */
    "degrees",                                       /* units           */
    -20.1, 20.1,                                     /* valid range     */
    0,0,                                             /* slope           */
    DFNT_FLOAT32,                                    /* HDF number type */
    1,                                               /* rank            */
    numScans, 1, 1,                                  /* dimension sizes */
    "Number of Scan Lines", NULL, NULL               /* dimension names */
    ) );

    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "slon",                                          /* short name      */
    "Scan start-pixel longitude",                    /* long name       */
    NULL,                                            /* standard name   */
    "degrees",                                       /* units           */
    -180.0, 180.0,                                   /* valid range     */
    0,0,                                             /* slope           */
    DFNT_FLOAT32,                                    /* HDF number type */
    1,                                               /* rank            */
    numScans, 1, 1,                                  /* dimension sizes */
    "Number of Scan Lines", NULL, NULL               /* dimension names */
    ) );

    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "clon",                                          /* short name      */
    "Scan center-pixel longitude",                   /* long name       */
    NULL,                                            /* standard name   */
    "degrees",                                       /* units           */
    -180.0, 180.0,                                   /* valid range     */
    0,0,                                             /* slope           */
    DFNT_FLOAT32,                                    /* HDF number type */
    1,                                               /* rank            */
    numScans, 1, 1,                                  /* dimension sizes */
    "Number of Scan Lines", NULL, NULL               /* dimension names */
    ) );

    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "elon",                                          /* short name      */
    "Scan end-pixel longitude",                      /* long name       */
    NULL,                                            /* standard name   */
    "degrees",                                       /* units           */
    -180.0, 180.0,                                   /* valid range     */
    0,0,                                             /* slope           */
    DFNT_FLOAT32,                                    /* HDF number type */
    1,                                               /* rank            */
    numScans, 1, 1,                                  /* dimension sizes */
    "Number of Scan Lines", NULL, NULL               /* dimension names */
    ) );

    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "slat",                                          /* short name      */
    "Scan start-pixel latitude",                     /* long name       */
    NULL,                                            /* standard name   */
    "degrees",                                       /* units           */
    -90.0, 90.0,                                     /* valid range     */
    0,0,                                             /* slope           */
    DFNT_FLOAT32,                                    /* HDF number type */
    1,                                               /* rank            */
    numScans, 1, 1,                                  /* dimension sizes */
    "Number of Scan Lines", NULL, NULL               /* dimension names */
    ) );

    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "clat",                                          /* short name      */
    "Scan center-pixel latitude",                    /* long name       */
    NULL,                                            /* standard name   */
    "degrees",                                       /* units           */
    -90.0, 90.0,                                     /* valid range     */
    0,0,                                             /* slope           */
    DFNT_FLOAT32,                                    /* HDF number type */
    1,                                               /* rank            */
    numScans, 1, 1,                                  /* dimension sizes */
    "Number of Scan Lines", NULL, NULL               /* dimension names */
    ) );

    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "elat",                                          /* short name      */
    "Scan end-pixel latitude",                       /* long name       */
    NULL,                                            /* standard name   */
    "degrees",                                       /* units           */
    -90.0, 90.0,                                     /* valid range     */
    0,0,                                             /* slope           */
    DFNT_FLOAT32,                                    /* HDF number type */
    1,                                               /* rank            */
    numScans, 1, 1,                                  /* dimension sizes */
    "Number of Scan Lines", NULL, NULL               /* dimension names */
    ) );

    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "csol_z",                                        /* short name      */
    "Scan center-pixel solar zenith angle",          /* long name       */
    NULL,                                            /* standard name   */
    "degrees",                                       /* units           */
    0.0, 180.0,                                      /* valid range     */
    0,0,                                             /* slope           */
    DFNT_FLOAT32,                                    /* HDF number type */
    1,                                               /* rank            */
    numScans, 1, 1,                                  /* dimension sizes */
    "Number of Scan Lines", NULL, NULL               /* dimension names */
    ) );

    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "longitude",                                     /* short name      */
    "Longitudes at control points",                  /* long name       */
    NULL,                                            /* standard name   */
    "degrees",                                       /* units           */
    -180.0, 180.0,                                   /* valid range     */
    0,0,                                             /* slope           */
    DFNT_FLOAT32,                                    /* HDF number type */
    2,                                               /* rank            */
    numScans, nctl, 1,                               /* dimension sizes */
    "Number of Scan Lines", "Number of Pixel Control Points", NULL
    ) );
 
    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "latitude",                                      /* short name      */
    "Latitudes at control points",                   /* long name       */
    NULL,                                            /* standard name   */
    "degrees",                                       /* units           */
    -90.0, 90.0,                                     /* valid range     */
    0,0,                                             /* slope           */
    DFNT_FLOAT32,                                    /* HDF number type */
    2,                                               /* rank            */
    numScans, nctl, 1,                               /* dimension sizes */
    "Number of Scan Lines", "Number of Pixel Control Points", NULL
    ) );

    /*                                                                  */
    /* Create the geophysical SDSes                                     */
    /* ---------------------------------------------------------------- */
    /*                                                                  */
    for (i=0; i<tot_prod; i++) {

        if ((p = get_l2prod_index(l2file->l2_prod_names[i],l2file->sensorID,
               numBands+numBandsIR,numPixels,numScans,Lambda)) == NULL) {
            fHY1C_out(stderr,
            "-E- %s line %d: product index failure.\n",
            __FILE__,__LINE__);
            return(1);
        };

        PTB( CreateSDS(
        sd_id,
        l2file->l2_prod_names[i],                    /* short name      */
        p->title,                                    /* long name       */
        p->standard_name,                            /* standard name   */
        p->units,                                    /* units           */
        p->min, p->max,                              /* valid range     */
        p->slope,p->offset,                          /* slope, offset   */
        p->datatype,                                 /* HDF number type */
        p->rank,                                     /* number of dims  */
        p->dim[0],p->dim[1],p->dim[2],               /* dimension sizes */
        p->dimname[0],p->dimname[1],p->dimname[2]    /* dimension names */
        ) );


        /*                                                              */
        /* Add solar irradiance meta-data if this is the nLw or Rrs     */
        /*                                                              */
        switch (p->cat_ix) {
          case CAT_nLw:
          case CAT_Rrs:
            sds_id = SDselect(sd_id, SDnametoindex(sd_id,l2file->l2_prod_names[i]));
            PTB( sd_setattr( sds_id,"solar_irradiance", DFNT_FLOAT32, 1, &Fonom[p->prod_ix]));
            SDendaccess(sds_id);
        }

        /*                                                              */
        /* Add flag-name attributes if this is the l2_flags product     */
        /*                                                              */
        if (strcmp(l2file->l2_prod_names[i],"l2_flags") == 0) {
            sds_id = SDselect(sd_id, SDnametoindex(sd_id, "l2_flags"));
            for (k=0; k<NFLAGS; k++) {
                PTB( sd_setattr( sds_id,l2_flag_sname[k], DFNT_CHAR,
                                 strlen(l2_flag_lname[k])+1,l2_flag_lname[k]) );
            }
            SDendaccess(sds_id);

        } else {

            /*                                                            */
            /* Add bad value attributes                                   */
            /*                                                            */
            sds_id = SDselect(sd_id, SDnametoindex(sd_id,l2file->l2_prod_names[i]));
            switch (p->datatype) {
              case DFNT_UINT8:
                pbuf = (VOIDP) float2uint8((VOIDP)&bad_float,0,1,1,p->slope,p->offset);
                sd_setattr( sds_id,"bad_value_scaled", DFNT_INT8,    1, pbuf);
                pbuf = (VOIDP) unscale_sds(pbuf,p,0,1,1);
                sd_setattr( sds_id,"bad_value_unscaled", DFNT_FLOAT32, 1, pbuf);         
                break;
              case DFNT_INT16:
	        pbuf = (VOIDP) float2int16((VOIDP)&bad_float,0,1,1,p->slope,p->offset);
                sd_setattr( sds_id,"bad_value_scaled", DFNT_INT16,   1, pbuf);         
                pbuf = (VOIDP) unscale_sds(pbuf,p,0,1,1);
                sd_setattr( sds_id,"bad_value_unscaled", DFNT_FLOAT32, 1, pbuf);         
                break;
              case DFNT_INT32:
                break;
              case DFNT_FLOAT32:
                pbuf = (VOIDP) float2float((VOIDP)&bad_float,0,1,1,p->slope,p->offset);
                sd_setattr( sds_id,"bad_value_scaled", DFNT_FLOAT32, 1, pbuf);         
                pbuf = (VOIDP) unscale_sds(pbuf,p,0,1,1);
                sd_setattr( sds_id,"bad_value_unscaled", DFNT_FLOAT32, 1, pbuf);         
                break;
              default:
                break;
            }
	}
    }


    /*                                                                  */
    /* Create the wavelength and control points array SDSes             */
    /* ---------------------------------------------------------------- */
    /*                                                                  */
    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "wavelength",                                    /* short name      */
    "Wavelengths",                                   /* long name       */
    NULL,                                            /* standard name   */
    "nm",                                            /* units           */
    0,0,                                             /* range           */
    0,0,                                             /* slope           */
    DFNT_INT32,                                      /* HDF number type */
    1,                                               /* rank            */
    numBands+numBandsIR, 0, 0,                       /* dimension sizes */
    "total band number", NULL, NULL                  /* dimension names */
    ) );

    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "cntl_pt_cols",                                  /* short name      */
    "Pixel control points",                          /* long name       */
    NULL,                                            /* standard name   */
    "none",                                          /* units           */
    0,0,                                             /* range           */
    0,0,                                             /* slope           */
    DFNT_INT32,                                      /* HDF number type */
    1,                                               /* rank            */
    nctl, 1, 1,                                      /* dimension sizes */
    "Number of Pixel Control Points",NULL, NULL      /* dimension names */
    ) );

    PTB( CreateSDS(
    sd_id,                                           /* file id         */
    "cntl_pt_rows",                                  /* short name      */
    "Scan control points",                           /* long name       */
    NULL,                                            /* standard name   */
    "none",                                          /* units           */
    0,0,                                             /* range           */
    0,0,                                             /* slope           */
    DFNT_INT32,                                      /* HDF number type */
    1,                                               /* rank            */
    numScans, 1, 1,                                  /* dimension sizes */
    "Number of Scan Lines",NULL, NULL                /* dimension names */
    ) );

    if (numBands > 0) {

        PTB( CreateSDS(
        sd_id,                                           /* file id         */
        "vcal_gain",                                     /* short name      */
        "Vicarious Calibration Gain",                    /* long name       */
        NULL,                                            /* standard name   */
        "dimensionless",                                 /* units           */
        0,0,                                             /* range           */
        0,0,                                             /* slope           */
        DFNT_FLOAT,                                      /* HDF number type */
        1,                                               /* rank            */
        numBands, 1, 1,                                  /* dimension sizes */
        "band number",NULL, NULL                         /* dimension names */
        ) );

    
        PTB( CreateSDS(
        sd_id,                                           /* file id         */
        "vcal_offset",                                   /* short name      */
        "Vicarious Calibration Offset",                  /* long name       */
        NULL,                                            /* standard name   */
        "mW cm^-2 um^-1 sr^-1",                          /* units           */
        0,0,                                             /* range           */
        0,0,                                             /* slope           */
        DFNT_FLOAT,                                      /* HDF number type */
        1,                                               /* rank            */
        numBands, 1, 1,                                  /* dimension sizes */
        "band number",NULL, NULL                         /* dimension names */
        ) );

    
        PTB( CreateSDS(
        sd_id,                                           /* file id         */
        "F0",                                            /* short name      */
        "Mean Solar Flux",                               /* long name       */
        NULL,                                            /* standard name   */
        "mW cm^-2 um^-1",                                /* units           */
        0,0,                                             /* range           */
        0,0,                                             /* slope           */
        DFNT_FLOAT,                                      /* HDF number type */
        1,                                               /* rank            */
        numBands, 1, 1,                                  /* dimension sizes */
        "band number",NULL, NULL                         /* dimension names */
        ) );

    
        PTB( CreateSDS(
        sd_id,                                           /* file id         */
        "k_oz",                                          /* short name      */
        "Ozone Absorption Coefficient",                  /* long name       */
        NULL,                                            /* standard name   */
        "cm^-1",                                         /* units           */
        0,0,                                             /* range           */
        0,0,                                             /* slope           */
        DFNT_FLOAT,                                      /* HDF number type */
        1,                                               /* rank            */
        numBands, 1, 1,                                  /* dimension sizes */
        "band number",NULL, NULL                         /* dimension names */
        ) );

    
        PTB( CreateSDS(
        sd_id,                                           /* file id         */
        "Tau_r",                                         /* short name      */
        "Rayleigh Optical Thickness",                    /* long name       */
        NULL,                                            /* standard name   */
        "dimensionless",                                 /* units           */
        0,0,                                             /* range           */
        0,0,                                             /* slope           */
        DFNT_FLOAT,                                      /* HDF number type */
        1,                                               /* rank            */
        numBands, 1, 1,                                  /* dimension sizes */
        "band number",NULL, NULL                         /* dimension names */
        ) );

    }

    
    /*                                                                  */
    /* Write out some global attributes                                 */
    /* ---------------------------------------------------------------- */
    /*                                                                  */
    sprintf(title,"%s Level-2 Data",sensorName[l2file->sensorID]);
    strcpy(soft_id, VERSION);

#ifdef SEADAS
    get_swf_def_meta1(seadas_vs, data_center, station_name, &station_lat,
                      &station_lon, station_code);
    strcpy(prog_version, PROGRAM);
    strcat(prog_version, " ");
    strcat(prog_version, VERSION);
    gen_soft_id(seadas_vs, prog_version, soft_id);
#endif

    PTB( SetChrGA (sd_id, "Title", title)                               );
    PTB( SetChrGA (sd_id, "Sensor Name" , sensorName[l2file->sensorID]) );
    PTB( SetChrGA (sd_id, "Product Name", basename(l2file->name))       );
    PTB( SetChrGA (sd_id, "Software Name", PROGRAM)                     );
    PTB( SetChrGA (sd_id, "Software Version", soft_id)                  );
    PTB( SetChrGA (sd_id, "Processing Version", l2file->input->pversion));

#ifdef SEADAS
    PTB( SetChrGA (sd_id, "Data Center", data_center)                   );
    PTB( SetChrGA (sd_id, "Station Name", station_name)                 );
    PTB( SetF32GA (sd_id, "Station Latitude", station_lat)              );
    PTB( SetF32GA (sd_id, "Station Longitude", station_lon)             );
#endif

    if (l2file->sensorID != SEAWIFS) {
     PTB( SetF32GA (sd_id, "Orbit Node Longitude", l2file->orbit_node_lon)  );
     PTB( SetI32GA (sd_id, "Orbit Number", l2file->orbit_number)  );
     PTB( SetChrGA (sd_id, "Node Crossing Time", l2file->node_crossing_time)  );
    }
     

    PTB( SetChrGA (sd_id, "Processing Time", ydhmsf(now(),'G'))         );
    PTB( SetChrGA (sd_id, "Processing Control", l2file->pro_control)    );
    PTB( SetChrGA (sd_id, "Input Parameters", l2file->input_parms)      );
    PTB( SetChrGA (sd_id, "Input Files", l2file->input_files)           );
    PTB( SetChrGA (sd_id, "Calibration Data", l2file->calfile)          );
    PTB( SetChrGA (sd_id, "Mask Names", l2file->mask_names)             );
    PTB( SetI32GA (sd_id, "Number of Bands", numBands)                  );
    PTB( SetI32GA (sd_id, "Number of Scan Lines", numScans)             );
    PTB( SetI32GA (sd_id, "Pixels per Scan Line", numPixels)            );
    PTB( SetI32GA (sd_id, "Scene Center Scan Line", cscan)              );
    PTB( SetI32GA (sd_id, "Number of Scan Control Points",numScans)     );
    PTB( SetI32GA (sd_id, "Number of Pixel Control Points",nctl)        );

    /*                                                                  */
    /* Write out some global SDSs                                       */
    /* ---------------------------------------------------------------- */
    /*                                                                  */
    PTB( sd_writedata(sd_id, "wavelength",   Lambda, 0, 0, 0, numBands+numBandsIR, 0,0) );

    if (numBands > 0) {
        PTB( sd_writedata(sd_id, "vcal_gain",    Gain,   0, 0, 0, numBands, 0,0) );
        PTB( sd_writedata(sd_id, "vcal_offset",  Offset, 0, 0, 0, numBands, 0,0) );
        PTB( sd_writedata(sd_id, "F0",           Fobar,  0, 0, 0, numBands, 0,0) );
        PTB( sd_writedata(sd_id, "k_oz",         k_oz,   0, 0, 0, numBands, 0,0) );
        PTB( sd_writedata(sd_id, "Tau_r",        Tau_r,  0, 0, 0, numBands, 0,0) );
        PTB( sd_writedata(sd_id, "cntl_pt_cols", ictl,   0, 0, 0, nctl,     0,0) );
        PTB( sd_writedata(sd_id, "cntl_pt_rows", jctl,   0, 0, 0, numScans, 0,0) );
    }

    return(LIFE_IS_GOOD);
}


/*-------- -------------------------------------------------------- */
/* Update scan-line SDSs for the specified scan.                    */
/* ---------------------------------------------------------------- */

int writel2_hdf( filehandle *l2file, int32_t recnum, l2str *l2rec) 
{
    static float32 *lon  = NULL;
    static float32 *lat  = NULL;
    static int32   *buf  = NULL;
    static uint8   *cbuf = NULL;
    static float32 *fbuf = NULL; 
    static float    fsol = -1.0;
    VOIDP   pbuf;
    int32   sd_id    = l2file->sd_id;
    int32   i,j,ip;
    l2prodstr *p;
    uint32_t mask;

    if (recnum >= numScans) {
       fHY1C_out(stderr,"-W- %s line %d: ", __FILE__,__LINE__);
       fHY1C_out(stderr,"attempt to write rec %d of %d\n",recnum,numScans);
       return(1);
    } 

    if (recnum >= cscan && fsol < 0.0) {
        fsol = l2rec->fsol;
    }

    /* allocate buffer space */
    if (buf == NULL) {
        if ( (buf = calloc(numPixels,sizeof(int32))) == NULL) {
            fHY1C_out(stderr,
            "-E- %s line %d: Unable to allocate buffer space.\n",
            __FILE__,__LINE__);
            exit(1);
        }
        cbuf = (uint8   *) buf;
        fbuf = (float32 *) buf;
    }
    if (lon == NULL) {
        if ( (lon = calloc(numPixels,sizeof(float32))) == NULL) {
            fHY1C_out(stderr,
            "-E- %s line %d: Unable to allocate buffer space.\n",
            __FILE__,__LINE__);
            exit(1);
        }
    }
    if (lat == NULL) {
        if ( (lat = calloc(numPixels,sizeof(float32))) == NULL) {
            fHY1C_out(stderr,
            "-E- %s line %d: Unable to allocate buffer space.\n",
            __FILE__,__LINE__);
            exit(1);
        }
    }


    /* Set arrays that depend on the number of control points */
    for (i = 0; i < nctl; i++) {
        lon [i] = l2rec->lon  [ictl[i]-1];
        lat [i] = l2rec->lat  [ictl[i]-1];
    }

    /* Write the scan-line data */
    PTB( sd_writedata(sd_id, "year", l2rec->year, recnum, 0,0,1,1,1) );
    PTB( sd_writedata(sd_id, "day" , l2rec->day , recnum, 0,0,1,1,1) );
    PTB( sd_writedata(sd_id, "msec", l2rec->msec, recnum, 0,0,1,1,1) );
    PTB( sd_writedata(sd_id, "tilt",&(l2rec->tilt), recnum, 0, 0, 1, 1, 1) );
    PTB( sd_writedata(sd_id, "slon",&(l2rec->lon[spix]),recnum, 0,0,1,1,1) );
    PTB( sd_writedata(sd_id, "clon",&(l2rec->lon[cpix]),recnum, 0,0,1,1,1) );
    PTB( sd_writedata(sd_id, "elon",&(l2rec->lon[epix]),recnum, 0,0,1,1,1) );
    PTB( sd_writedata(sd_id, "slat",&(l2rec->lat[spix]),recnum, 0,0,1,1,1) );
    PTB( sd_writedata(sd_id, "clat",&(l2rec->lat[cpix]),recnum, 0,0,1,1,1) );
    PTB( sd_writedata(sd_id, "elat",&(l2rec->lat[epix]),recnum, 0,0,1,1,1) );
    PTB( sd_writedata(sd_id, "csol_z",&(l2rec->solz[cpix]),recnum, 0,0,1,1,1) );
    PTB( sd_writedata(sd_id, "longitude", lon, recnum, 0,0,1,nctl,1) );
    PTB( sd_writedata(sd_id, "latitude" , lat, recnum, 0,0,1,nctl,1) );


    /*                                                                  */
    /* Write the geophysical data                                       */
    /* ---------------------------------------------------------------- */
    /*                                                                  */
    for (i=0; i<l2file->tot_prod; i++) {

        // get the product index record

        if ((p = get_l2prod_index(l2file->l2_prod_names[i],l2file->sensorID,
                 numBands+numBandsIR,numPixels,numScans,Lambda)) == NULL) {
            fHY1C_out(stderr,
            "-E- %s line %d: product index failure.\n",
            __FILE__,__LINE__);
            return(1);
        };

        // extract the product and scale (if needed)

        if (p->slope == 1.0 && p->offset == 0.0)
            pbuf = prodgen(p,l2rec);
        else
            pbuf = scale_sds( (float *) prodgen(p,l2rec), p, 1);

        // write to L2 file

        if (p->rank == 2) {
            PTB( sd_writedata(sd_id,l2file->l2_prod_names[i],pbuf,
                              recnum,0,0,1,numPixels,1) );
	} else {
            PTB( sd_writedata(sd_id,l2file->l2_prod_names[i],pbuf,
                              recnum,0,0,1,1,1) );
	}       
    }

    /* Update flag counter */
    for (ip=0; ip<l2rec->npix; ip++) {
        mask = 1L;
        for (i=0; i<NFLAGS; i++) {
            l2file->flag_cnt[i] += ((l2rec->flags[ip] & mask) > 0);
            mask *= 2L;
	}
    }

    /* Write global attributes */
    if (recnum == (numScans- 1)) {

        float flag_perc[NFLAGS];

        scene_meta_write(sd_id);    

        PTB( SetF64GA(sd_id, "Earth-Sun Distance Correction",fsol));
    
        /* Report flag percentages */
        HY1C_out("\nPercentage of pixels flagged:\n");
        for (i=0; i<NFLAGS; i++) {
            flag_perc[i] = ((float) l2file->flag_cnt[i])/numScans/numPixels*100.0;
            HY1C_out("Flag #%2d: %16s %10d %8.4f\n",
                i+1,l2_flag_lname[i],l2file->flag_cnt[i],flag_perc[i]);
            if (fp_meta != NULL) 
                fHY1C_out(fp_meta,"Flag #%2d: %16s %10d %8.4f\n",
                    i+1,l2_flag_lname[i],l2file->flag_cnt[i],flag_perc[i]);
	}
        PTB(sd_setattr(sd_id,"Flag Percentages",
            DFNT_FLOAT32,NFLAGS,(VOIDP)flag_perc));

    }   

    return(LIFE_IS_GOOD);
}

/* -------------------------------------------------------- */
/* Finish access for the current file.                      */
/* -------------------------------------------------------- */
int closel2_hdf(filehandle *l2file)
{
    PTB( MakeVgroups(l2file) );

    if (SDend(l2file->sd_id)) {
       fHY1C_out(stderr,"-E- %s line %d: SDend(%d) failed for file, %s.\n",
       __FILE__,__LINE__,l2file->sd_id,l2file->name);
       return(HDF_FUNCTION_ERROR);
    }

    if (fp_meta != NULL)
        fclose(fp_meta);

    return(LIFE_IS_GOOD);
}
